//============================================================================
// Name        : FabVoronoi.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <iterator>
#include <iomanip>
#include "Topol.h"
#include "FAtoms.h"
#include "HeaderTrj.h"
#include "ResidueCM.h"
#include "FVoronoi.h"
#include "FComms.h"
#include "Help.hpp"
#include "FstreamC.h"
#include "FstreamF.h"
#include "VoronoiPrint.h"
#include "trjInput.h"

static int pGroup=-1; // From 0 to ngrps, ngrps no. of selected groups
static bool bPrintVols=true;
static bool bPrintAreas=FALSE;


using namespace Topol_NS;
using namespace std;
int main(int argc, char ** argv)
{
	trj::trjInput MyInput(argc,argv);
    Parallel::comm=new Parallel::MyMPI(argc,argv);
    bool bHyd=true;
	bool bDel=true;
	bool bWriteFab=false;
	bool bIsrd=false;
	bool bTestVol=false;
	int narg=1;
	string filein,fileout,filepdb;
	string scratchdir;
	ifstream fpdb;
	Fstream * finx;
	ofstream fout;
	string comm0(argv[0]);
	size_t mypos=comm0.find_last_of("/")+1;
	size_t length=comm0.size()-mypos;
	string command=comm0.substr(mypos,length);
	string errmsg;
	string Usage="Usage:\t"+ command + "\n";
	vector<string> use=MyInput.getUsage();
	for(unsigned int n=0;n<use.size();n++)
		Usage+=use[n];
	Usage+="\n\t Default values in square brackets []\n";
	try{
		if(argc == 1) throw Usage;
		if(int m=MyInput.bTest().size()) {
			errmsg=" Command(s) not found: ";
			for(unsigned int n=0;n<m;n++)
				errmsg+=MyInput.bTest()[n]+"  ";
			errmsg+="\n"+Usage;
			throw errmsg;
			}
		}
	catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}



	int nstart=0;
	int nend=-1;
	int nskip=1;
	// input read
	try{
		if(!MyInput["-dcd"].empty()) {
			if(MyInput["-dcd"].size() != 2) throw string(" filename expected for " + MyInput["-dcd"][0] + " option ");
			filein=MyInput["-dcd"][1];
			finx=new FstreamF(filein);
		}
		if(!MyInput["-xtc"].empty()) {
			if(MyInput["-xtc"].size() != 2) throw string(" filename expected for " + MyInput["-xtc"][0] + " option ");
			filein=MyInput["-xtc"][1];
			finx=new FstreamC(filein);
		}
		if(!MyInput["-pdb"].empty()) {
			if(MyInput["-pdb"].size() != 2) throw string(" filename expected for " + MyInput["-pdb"][0] + " option ");
			filepdb=MyInput["-pdb"][1];
			fpdb.open(filepdb.c_str(),ios::in);
		}
		if(!MyInput["-o"].empty()) {
			if(MyInput["-o"].size() != 2) throw string(" filename expected for " + MyInput["-o"][0] + " option ");
			fileout=MyInput["-o"][1];
		}
		if(!MyInput["-b"].empty()) {
			if(MyInput["-b"].size() != 2) throw string(" Number of first frame needed for " + MyInput["-b"][0] + " option ");
			stringstream(MyInput["-b"][1])>> nstart;
		}
		if(!MyInput["-e"].empty()) {
			if(MyInput["-e"].size() != 2) throw string(" Number of last frame needed for " + MyInput["-e"][0] + " option ");
			stringstream(MyInput["-e"][1])>> nend;
		}
		if(!MyInput["-skip"].empty()) {
			if(MyInput["-skip"].size() != 2) throw string(" Number of skipped frames needed for " + MyInput["-skip"][0] + " option ");
			stringstream(MyInput["-skip"][1])>> nskip;
		}
		if(!MyInput["-pvol"].empty()) {
			if(MyInput["-pvol"].size() != 1) throw string(" No argument for " + MyInput["-pvol"][0] + " option ");
				bPrintVols=true;
		}
		if(!MyInput["-pgroup"].empty()) {
			if(MyInput["-pgroup"].size() != 2) throw string(" Group number needed for " + MyInput["-pgroup"][0] + " option ");
			stringstream(MyInput["-pgroup"][1])>> pGroup;
		}
		if(!MyInput["-parea"].empty()) {
			if(MyInput["-parea"].size() != 1) throw string(" No argument for " + MyInput["-parea"][0] + " option ");
				bPrintAreas=true;
		}
		if(!MyInput["-nohyd"].empty()) {
			if(MyInput["-nohyd"].size() != 1) throw string(" No argument for " + MyInput["-nohyd"][0] + " option ");
				bHyd=false;
		}
		if(!MyInput["-hyd"].empty()) {
			if(MyInput["-hyd"].size() != 1) throw string(" No argument for " + MyInput["-hyd"][0] + " option ");
				bHyd=true;
		}
		if(!MyInput["-nodel"].empty()) {
			if(MyInput["-nodel"].size() != 1) throw string(" No argument for " + MyInput["-nodel"][0] + " option ");
				bDel=false;
		}
		if(!MyInput["-del"].empty()) {
			if(MyInput["-del"].size() != 1) throw string(" No argument for " + MyInput["-del"][0] + " option ");
				bDel=false;
		}
		if(!MyInput["-rd"].empty()) {
			if(MyInput["-rd"].size() != 1) throw string(" No argument for " + MyInput["-rd"][0] + " option ");
				bIsrd=true;
		}
		if(!MyInput["-nord"].empty()) {
			if(MyInput["-nord"].size() != 1) throw string(" No argument for " + MyInput["-nord"][0] + " option ");
				bIsrd=false;
		}
		if(!MyInput["-fab"].empty()) {
			if(MyInput["-fab"].size() != 1) throw string(" No argument for " + MyInput["-fab"][0] + " option ");
				bWriteFab=true;
		}
		if(!MyInput["-nofab"].empty()) {
			if(MyInput["-nofab"].size() != 1) throw string(" No argument for " + MyInput["-nofab"][0] + " option ");
			bWriteFab=false;
		}
		if(!MyInput["-test"].empty()) {
			if(MyInput["-test"].size() != 1) throw string(" No argument for " + MyInput["-test"][0] + " option ");
				bTestVol=true;
		}
		if(!MyInput["-help"].empty()) {
			if(MyInput["-help"].size() != 1) throw string(" No argument for " + MyInput["-help"][0] + " option ");
				throw Help;
		}
	}catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}
	VoronoiPrint::bPrintAreas=bPrintAreas;
    VoronoiPrint::bPrintVols=bPrintVols;
    VoronoiPrint::pGroup=pGroup;

	Fstream & fin=*finx;

// read pdb file to construct topology
	vector<string> data;
	for(string str;getline(fpdb,str);){
		data.push_back(str);
	}

	Topol MyTop(data,bIsrd);

	int natoms=MyTop.Size();
	int nresid=MyTop.ResSize();
	FAtoms atm;
	atm.setDim(natoms);

// Map atoms to residues to use in Voronoi. Eliminate hydrogens if needed
	vector<int> * resmap=new vector<int> [nresid];
    string * ResLabels=new string [nresid];
    string label;

    for(int i=0;i<nresid;i++){
    	ResLabels[i]=MyTop.ResInfo(i);
    }
    for(int i=0;i<natoms;i++){
    	int j=MyTop.ResIndex(i);

    	if(bHyd) resmap[j].push_back(i);
    	else {
    		label=MyTop.AtomInfo(i);
    		if(label.compare(0,1,"H",0,1) != 0) resmap[j].push_back(i);
    	}
    }
    ResidueCM MyCM;
    MyCM(natoms,nresid,resmap,ResLabels);

    unsigned int totaltypes=MyTop.getCidxSize();
    AtomIndex * cidx=new AtomIndex [totaltypes];
    for(unsigned int i=0;i<MyTop.getCidxSize();i++){
		const vector<int> & ty=MyTop.getCidx(i);
		cidx[i].set(ty.size(),&ty[0]);
	}

	MyCM.ResidueSelect(totaltypes,cidx);
// Setup Voronoi

	FVoronoi * vor=new FVoronoi(natoms,MyTop.getAtomName(),bHyd);
	vor->setTypes(MyTop.getmaxt(),cidx,MyTop.getTypeNames());
	Voronoi * RegVor=dynamic_cast<Voronoi *> (vor);

// Read header of dcd file
	fin.seekg(0,"end");
	ios::streampos len=fin.tellg();
	fin.seekg(0,"beg");
	HeaderTrj header;
	fin>>header;

	try{if(!header.check(MyTop.Size())) throw string("Number of atoms in the pdb and trajectory files does not match!");}
	catch(const string & s){cout << s<<endl;exit(-1);}

	atm.setrd(MyTop.getrd());

// Read trajectory and compute Voronoi data

    Parallel::comm->Barrier();


	Parallel::FComms Comms(fout,fileout,nstart,nend,header.getNFR(),nskip);
	nstart=Comms.getStart();
	nend=Comms.getEnd();
	int ntime=0;

	while(fin.tellg() < len){
		if(ntime < nstart) {fin+=atm;ntime++;continue;}
		if(nend != -1 && ntime > nend) break;
		fin >> atm;
		cout << fixed << setw(5) << "----> Step No. " << ntime <<"\n";
		vor->Start(static_cast<float> (ntime),atm);
		vor->getData();
		if(bWriteFab)
			Comms.getStream() << *vor;
		else
			Comms.getStream() << *RegVor;
		if(bTestVol) vor->testVol();
		ntime++;
	};

	Comms.appendStreams();
	if(bDel) Comms.removeFiles();
	Parallel::comm->Barrier();
	Parallel::comm->~MyMPI();
	cout << "\nProgram completed: Output data written to " + fileout << "\n\n";
    return 0;
}   //end of main function
